#!/usr/bin/env python
############################################################################
#
# MODULE:       r.mcda.roughset
# AUTHOR:       Gianluca Massei - Antonio Boggia
# PURPOSE:      Generate a MCDA map from several criteria maps using Dominance Rough Set Approach - DRSA
#                       (DOMLEM algorithm proposed by  (S. Greco, B. Matarazzo, R. Slowinski)
# COPYRIGHT:  c) 2010 Gianluca Massei, Antonio Boggia  and the GRASS
#                       Development Team. This program is free software under the
#                       GNU General PublicLicense (>=v2). Read the file COPYING
#                       that comes with GRASS for details.
#
#############################################################################


import sys
import copy
import numpy as np
from time import time, ctime


def BuildFileISF(attributes, preferences, decision, outputMap, outputTxt):
    outputTxt=outputTxt+".isf"
    outf = file(outputTxt,"w")
    outf.write("**ATTRIBUTES\n")
    for i in range(len(attributes)):
        outf.write("+ %s: (continuous)\n" % attributes[i])
    outf.write("+ %s: [" % decision)
    value=[]
    value=grass.read_command("r.describe", flags = "1n", map = decision)
    v=value.split()

    for i in range(len(v)-1):
        outf.write("%s, " % str(v[i]))
    outf.write("%s]\n" % str(v[len(v)-1]))
    outf.write("decision: %s\n" % decision)

    outf.write("\n**PREFERENCES\n")
    for i in range(len(attributes)):
        if(preferences[i]==""):
            preferences[i]="none"
        outf.write("%s: %s\n" % (attributes[i], preferences[i]))
    outf.write("%s: gain\n" % decision)


    outf.write("\n**EXAMPLES\n")
    examples=[]
    MATRIX=[]
    for i in range(len(attributes)):
        grass.mapcalc("rast=if(isnull(${decision})!=0,${attribute},null())",
                         rast="rast",
                         decision=decision,
                         attribute=attributes[i])
        tmp=grass.read_command("r.stats", flags = "1n", nv="?", input = "rast")
        example=tmp.split()
        examples.append(example)

    tmp=grass.read_command("r.stats", flags = "1n", nv="?", input = decision)
    example=tmp.split()
    examples.append(example)

    MATRIX=map(list,zip(*examples))
    MATRIX=[r for r in MATRIX if not '?' in r] #remove all rows with almost one "?"
    MATRIX=[list(i) for i in set(tuple(j) for j in MATRIX)] #remove duplicate example

    for r in range(len(MATRIX)):
        for c in range(len(MATRIX[0])):
            outf.write("%s " % round(float(MATRIX[r][c]), 2))
        outf.write("\n")

    outf.write("**END")
    outf.close()
    return outputTxt



def collect_attributes (data):
    "Collects the values of header files isf, puts them in an array of dictionaries"
    header=[]
    attribute=dict()
    j=0
    start=(data.index(['**ATTRIBUTES'])+1)
    end=(data.index(['**PREFERENCES'])-1)
    for r in range(start, end):
        attribute={'name':data[r][1].strip('+:')}
        header.append(attribute)
    decision=data[end-1][1]
    end=(data.index(['**EXAMPLES']))

    start=(data.index(['**PREFERENCES'])+1)
    for r in header:
        r['preference']=data[start+j][1]
        j=j+1
    return header


def collect_examples (data):
    "Collect examples values and put them in a matrix (list of lists) "

    matrix=[]
    data=[r for r in data if not '?' in r] #filter objects with " ?"
#    data=[data.remove(r) for r in data if data.count(r)>1]
    start=(data.index(['**EXAMPLES'])+1)
    end=data.index(['**END'])
    for i in range(start, end):
        data[i]=(map(float, data[i]))
        matrix.append(data[i])
    i=1
    for r in matrix:
        r.insert(0, str(i))
        i=i+1
##    matrix=[list(i) for i in set(tuple(j) for j in matrix)] #remove duplicate example
    return matrix


def FileToInfoSystem(isf):
    "Read *.isf file and copy it's values in Infosystem dictionary"
    data=[]
    try:
        infile=open(isf,"r")
        rows=infile.readlines()
        for line in rows:
            line=(line.split())
            if (len(line)>0 ):
                data.append(line)
        infile.close()
        infosystem={'attributes':collect_attributes(data),'examples':collect_examples(data)}
    except TypeError:
        print "\n\n Computing error or input file %s is not readeable. Exiting gracefully" % isf
        sys.exit(0)

    return infosystem


def UnionOfClasses (infosystem):
    "Find upward and downward union for all classes and put it in a dictionary"
    DecisionClass=[]
    AllClasses=[]
    matrix=infosystem['examples']
    for r in matrix:
        DecisionClass.append(int(r[-1]))
    DecisionClass=list(set(DecisionClass))
    for c in range(len(DecisionClass)):
        tmplist=[r for r in matrix if int(r[-1])==DecisionClass[c]]
        AllClasses.append(tmplist)

    return AllClasses


def DownwardUnionsOfClasses (infosystem):
    "For each decision class, downward union corresponding to a decision class\
    is composed of this class and all worse classes (<=)"

    DownwardUnionClass=[]
    DecisionClass=[]
    matrix=infosystem['examples']
    for r in matrix:
        DecisionClass.append(int(r[-1]))
    DecisionClass=list(set(DecisionClass))
    for c in DecisionClass:
        tmplist=[r for r in matrix if int(r[-1])<=c]
        DownwardUnionClass.append(tmplist)
        #label=[row[0] for row in tmplist]
    return DownwardUnionClass


def UpwardUnionsOfClasses (infosystem):
    "For each decision class, upward union corresponding to a decision class \
    is composed of this class and all better classes.(>=)"

    UpwardUnionClass=[]
    DecisionClass=[]
    matrix=infosystem['examples']
    for r in matrix:
        DecisionClass.append(int(r[-1]))
    DecisionClass=list(set(DecisionClass))
    for c in DecisionClass:
        tmplist=[r for r in matrix if int(r[-1])>=c]
        UpwardUnionClass.append(tmplist)
        #label=[row[0] for row in tmplist]
    return UpwardUnionClass


###############################
def is_better (r1,r2, preference):
    "Check if r1 is better than r2"
    return all((( x >=y and p=='gain') or (x<=y and p=='cost')) for x,y, p in zip(r1,r2, preference) )


def is_worst (r1,r2, preference):
    "Check if r1 is worst than r2"
    return all((( x <=y and p=='gain') or (x>=y and p=='cost')) for x,y, p in zip(r1,r2, preference) )
 #################################


def DominatingSet (infosystem):
    "Find P-dominating set"
    matrix=infosystem['examples']
    preference=[s['preference'] for s in infosystem['attributes'] ]
    Dominating=[]
    for row in matrix:
        examples=[r  for r in matrix if  is_better(r[1:-1], row[1:-1], preference) ]
        Dominating.append({'object':row[0], 'dominance':[i[0] for i in examples], 'examples':examples})
##    for dom in Dominating:
##        print  dom['dominance'] ,' dominating ', dom['object']
    return Dominating

def DominatedSet (infosystem):
    "Find P-Dominated set"
    matrix=infosystem['examples']
    preference=[s['preference'] for s in infosystem['attributes'] ]
    Dominated=[]
    for row in matrix:
        examples=[r  for r in matrix if  is_worst(r[1:-1], row[1:-1], preference[:-1]) ]
        Dominated.append({'object':row[0], 'dominance':[i[0] for i in examples], 'examples':examples})
##    for dom in Dominated:
##        print  dom['dominance'] ,' is dominated by ', dom['object']
    return Dominated


def LowerApproximation (UnionClasses,  Dom):
    "Find Lower approximation and return a dictionaries list"
    c=1
    LowApprox=[]
    single=dict()
    for union in UnionClasses:
        tmp=[]
        UClass=set([row[0] for row in union] )
        for d in Dom:
            if (UClass.issuperset(set(d['dominance']))): #if Union class is a superse of dominating/dominated set, =>single Loer approx.
                tmp.append(d['object'])
        single={'class':c, 'objects':tmp} #dictionary for lower approximation  --
        LowApprox.append(single) #insert all Lower approximation in a list
        c+=1
    return LowApprox


def UpperApproximation (UnionClasses,  Dom):
    "Find Upper approximation and return a dictionaries list"
    c=1
    UppApprox=[]
    single=dict()
    for union in UnionClasses:
        UnClass=[row[0] for row in union]  #single union class
        s=[]
        for d in Dom:
           if len(set(d['dominance']) & set(UnClass)) >0:
               s.append(d['object'])
#               print set(s)
        single={'class':c,'objects':list(set(s))}
        UppApprox.append(single)
        c+=1

    return UppApprox


def Boundaries (UppApprox, LowApprox):
    "Find Boundaries like doubtful regions"
    Boundary=[]
    single=dict()

    for i in range(len(UppApprox)):
        single={'class':i, 'objects':list (set(UppApprox[i]['objects'])-set(LowApprox[i]['objects']) )}
        Boundary.append(single)

    return Boundary


def AccuracyOfApproximation(UppApprox, LowApprox):
    "Define the accuracy of approximation of Upward and downward approximation class"
    return len(LowApprox)/len(UppApprox)


def QualityOfQpproximation(DownwardBoundary,  infosystem):
    "Defines the quality of approximation of the partition Cl or, briefly, the quality of sorting"
    UnionBoundary=set()
    U=set([i[0] for i in infosystem['examples']])
    for b in DownwardBoundary:
        UnionBoundary=set(UnionBoundary) | set(b['objects'])
    return float(len(U-UnionBoundary)) / float(len(U))


def FindObjectCovered (rules, selected):
    "Find objects covered by a single rule and return\
    all related examples covered"
    obj=[]
    examples=[]

    for rule in rules:
        examples.append(rule['objectsCovered'])

    if len(examples)>0:
        examples = reduce(set.intersection,map(set,examples))  #functional approach: intersect all lists if example is not empty
        examples = list(set(examples) & set([r[0] for r in selected]))
    return examples #all examples covered from a single rule


def Evaluate (elem,rules,G,selected,infosystem):
    "Calcolate first and second evaluate index, according with original DOMLEM Algorithm"
    tmpRules=copy.deepcopy(rules)
    tmpElem=copy.deepcopy(elem)
    tmpRules.append(tmpElem)
    Object=[]
    Object=FindObjectCovered(tmpRules,selected)
    if(float(len(Object)))>0:
             firstEvaluate=float(len(set(G) & set(Object))) / float(len(Object))
             secondEvaluate=float(len(set(G) & set(Object)))
    else:
             firstEvaluate=0
             secondEvaluate=0

    return firstEvaluate,secondEvaluate



def FindBestCondition (best, elem, rules, selected, G, infosystem):
    "Choose the best condition"

    firstElem,secondElem=Evaluate(elem,rules,G,selected,infosystem)
    firstBest,secondBest=Evaluate(best,rules,G,selected,infosystem)

    if (firstElem>firstBest) or (firstElem==firstBest and secondElem>=secondBest):
        best=copy.deepcopy(elem)
    else:
        best=best

    return best


def CheckMinimalCondition (rules,B,matrix):
    "Check minimal elementary condition from each rule"
    obj_cov_by_rules=[]
    if len(rules)>1:
        for e in rules:
            check=copy.deepcopy(rules)
            check.remove(e)
            obj_cov_by_rules=FindObjectCovered (rules, matrix)
            if set(obj_cov_by_rules).issubset(set(B)):
                rules=check
    return rules

def Type_one_rule (c,  e,  preference,  matrix):
    elem={'criterion':c,'condition':e, 'sign':preference[c-1],'class':'', \
    'objectsCovered':[r[0] for r in matrix if (((r[c] >= e ) and (preference[c-1] == 'gain')) \
                                                                              or ((r[c] <= e ) and (preference[c-1] == 'cost' )))],'label':''}
    return elem

def Type_three_rule (c,  e,  preference,  matrix):
    elem={'criterion':c,'condition':e, 'sign':preference[c-1],'class':'', \
    'objectsCovered':[r[0] for r in matrix if (((r[c] <= e ) and (preference[c-1] == 'gain')) \
                                                        or ((r[c] >= e ) and (preference[c-1] == 'cost' )))],'label':''}
    return elem


def Find_rules (B, infosystem, type_rule):
    "Search rule from a family of lower approximation of upward unions \
    of decision classes"
    start=time()
    matrix=copy.deepcopy(infosystem['examples'])
    criteria_num=len(infosystem['attributes'])
    criteria=[r[1:-1] for r in matrix]
    preference=[s['preference'] for s in infosystem['attributes'] ] #extract preference label
    num_rules=0 #total rules number for each lower approximation
    G=copy.deepcopy(B)        #a set of objects from the given approximation
    E=[]        #a set  of rules covering set B (is a list of dictionary)
    all_obj_cov_by_rules=[] #all objects covered by all rules in E
    selected=copy.deepcopy(matrix) #storage reduct matrix by single elementary condition
    while (len(G)!=0 ):
        rules=[]     #starting comples (single rule built from elementary conditions  )
        S=copy.deepcopy(G)         #set of objects currently covered by rule
        control=0
        while (len(rules)==0 or set(obj_cov_by_rules).issubset(B)==False):
            obj_cov_by_rules=[] #set covered by rules
            best={'criterion':'','condition':'','sign':'','class':'','objectsCovered':'','label':'', 'type':''} #best candidate for elementary condition - start as empty
            for c in range(1, criteria_num):
                Cond=[r[c] for r in selected if  r[0] in S] #for each positive object from S create an elementary condition
                for e in Cond:
                    if type_rule=='one':
                        elem= Type_one_rule (c,  e,  preference,  matrix)
                    elif type_rule=='three':
                        elem= Type_three_rule (c,  e,  preference,  matrix)

                    else:
                        elem={'criterion':'','condition':'','sign':'','class':'','objectsCovered':'','label':'', 'type':''}

                    best=FindBestCondition(best, elem, rules, selected, G, infosystem)
            if best not in rules:
                rules.append(best)   #add the best condition to the complex

            for r in rules:
                obj_cov_by_rules.append(r['objectsCovered'])
            obj_cov_by_rules=list((reduce(set.intersection,map(set,obj_cov_by_rules)))) #reduce():Apply function of two arguments cumulatively to the items of iterable, from left to right, so as to reduce the iterable to a single value.

            S=list(set(S) & set(best['objectsCovered'] ))
            control+=1

#        rules=CheckMinimalCondition (rules,B,matrix)

        if rules not in E:
            E.append(rules) #add the induced rule
            num_rules+=1
        all_obj_cov_by_rules=list(set(all_obj_cov_by_rules) | set(obj_cov_by_rules))

        G=list(set(B)-set(all_obj_cov_by_rules)) #remove example coverred by all finded rule -- this operation is a set difference
        selected=[o  for o in selected if  not  o[0]  in all_obj_cov_by_rules] #reduct matrix, remove object coverred by all finded rule
        num_rules+=1

    return E



def Domlem(Lu,Ld, infosystem):
    "DOMLEM algoritm \
    (An algorithm for induction of decision rules consistent with the dominance\
    principle - Greco S., Matarazzo, B., Slowinski R., Stefanowski J.)"
    attributes=infosystem['attributes']

    RULES=[]
## *** AT LEAST {>= Class} - Type 1 rules *** "
    for a in Lu[1:]:
        B=a['objects']
        E=Find_rules (B, infosystem, 'one')
        for e in E:
            for i in e:
                i['class']=a['class']
                i['label']=attributes[i['criterion']-1]['name']
                i['type']='at_last'
                if (attributes[i['criterion']-1]['preference']=='gain'):
                    i['sign']='>='
                else:
                    i['sign']='<='
            RULES.append(e)


##  *** AT MOST {<= Class} - Type 3 rules ***"
    for b in Ld[:-1]:
        B=b['objects']
        E=Find_rules(B, infosystem, 'three')
        for e in E:
            for i in e:
                i['class']=b['class']
                i['label']=attributes[i['criterion']-1]['name']
                i['type']='at_most'
                if (attributes[i['criterion']-1]['preference']=='gain'):
                    i['sign']='<='
                else:
                    i['sign']='>='
            RULES.append(e)

    return RULES


def Print_rules(RULES, outputTxt):
    "Print rls output file"

    i=1
    outfile=open(outputTxt+".rls","w")
    outfile.write('[RULES]\n')

    for R in RULES:
        outfile.write("%d: " % i, )
        #"&".join(e)
        for e in R:
                outfile.write("( %s %s %.3f )" % (e['label'], e['sign'],e['condition'] ))
        outfile.write("=> ( class %s , %s )\n" % ( e['type'], e['class'] ))
        i+=1
    outfile.close()
    return 0

def Parser_mapcalc(RULES, outputMap):
    "Parser to build a formula to be included  in mapcalc command"
    i=1
    category=[]
    maps=[]
    stringa=[]
    for R in RULES:
        formula="if("
        for e in R[:-1]: #build a mapcalc formula
            formula+= "(%s %s %.4f ) && " % (e['label'],  e['sign'] ,  e['condition'] )
        formula+= "(%s %s %.4f ),%d,null())" % (R[-1]['label'],R[-1]['sign'], R[-1]['condition'] ,i )
        mappa="%d.%s_%d" % ( i, R[0]['type'], R[0]['class'] ) #build map name for mapcalc output
        category.append({'id':i, 'type': R[0]['type'], 'class':R[0]['class']}) #extract category name
        maps.append(mappa) #extract maps name
        grass.mapcalc(mappa +"=" +formula)
        i+=1
    maps=",".join(maps)
    grass.run_command("r.patch", input=maps,  flags='-o', output=outputMap)
#    grass.run_command("g.remove",  rast=maps)
    rls = file("rules","w")
    for i in category:
        rls.write("%d:%s %s\n" % (i['id'],  i['type'], str(i['class'])))
    rls.close
    grass.run_command("r.category", map='classify', rules="rules"  )
    return 0

def main():
    "main function for DOMLEM algorithm"

    start=time()
    attributes = options['criteria'].split(',')
    preferences=options['preferences'].split(',')
    decision=options['decision']
    outputMap= options['outputMap']
    outputTxt= options['outputTxt']
    out=BuildFileISF(attributes, preferences, decision, outputMap, outputTxt)

    infosystem=FileToInfoSystem(out)

    UnionOfClasses(infosystem)
    DownwardUnionClass=DownwardUnionsOfClasses(infosystem)
    UpwardUnionClass=UpwardUnionsOfClasses(infosystem)
    Dominating=DominatingSet(infosystem)
    Dominated=DominatedSet(infosystem)
##    upward union class
    print "elaborate upward union"
    Lu=LowerApproximation(UpwardUnionClass, Dominating) #lower approximation of upward union for type 1 rules
    Uu=UpperApproximation(UpwardUnionClass,Dominated ) #upper approximation of upward union
    UpwardBoundary=Boundaries(Uu, Lu)
##    downward union class
    print "elaborate downward union"
    Ld=LowerApproximation(DownwardUnionClass, Dominated) # lower approximation of  downward union for type 3 rules
    Ud=UpperApproximation(DownwardUnionClass,Dominating ) # upper approximation of  downward union
    DownwardBoundary=Boundaries(Ud, Ld)

    QualityOfQpproximation(DownwardBoundary,  infosystem)
    print "RULES extraction"
    RULES=Domlem(Lu,Ld, infosystem)

    Parser_mapcalc(RULES, outputMap)
    Print_rules(RULES, outputTxt)
    end=time()
    print "Time computing-> %.4f s" % (end-start)

    return 0

###########execute the script##########################
if __name__=='__main__':
    main()
#######################################################

